
make_regs_bw_plus <- function(bw){

  df <- 
    individual_level %>% 
    mutate(forcing = (forcing + 10000),
           opium_legal = case_when(forcing > 0 ~ 1,
                                   TRUE ~ 0)) %>%
    filter(abs(forcing) < bw)
  
  mod1 = lm_robust(outcome1 ~ opium_legal + forcing + opium_legal*forcing, data = df, clusters = mkid14) %>% tidy()
  mod2 = lm_robust(outcome2 ~ opium_legal + forcing + opium_legal*forcing, data = df, clusters = mkid14) %>% tidy()
  mod3 = lm_robust(outcome3 ~ opium_legal + forcing + opium_legal*forcing, data = df, clusters = mkid14) %>% tidy()
  mod4 = lm_robust(outcome4 ~ opium_legal + forcing + opium_legal*forcing, data = df, clusters = mkid14) %>% tidy()
  
  bind_rows(mod1, mod2, mod3, mod4) %>%
    mutate(bw = bw,
           placebo = 10000) %>%
    filter(term == "opium_legal") %>%
    return()
}

make_regs_bw_minus <- function(bw){
  
  df <- 
    individual_level %>% 
    mutate(forcing = (forcing - 10000),
           opium_legal = case_when(forcing > 0 ~ 1,
                                   TRUE ~ 0)) %>%
    filter(abs(forcing) < bw)
  
  mod1 = lm_robust(outcome1 ~ opium_legal + forcing + opium_legal*forcing, data = df, clusters = mkid14) %>% tidy()
  mod2 = lm_robust(outcome2 ~ opium_legal + forcing + opium_legal*forcing, data = df, clusters = mkid14) %>% tidy()
  mod3 = lm_robust(outcome3 ~ opium_legal + forcing + opium_legal*forcing, data = df, clusters = mkid14) %>% tidy()
  mod4 = lm_robust(outcome4 ~ opium_legal + forcing + opium_legal*forcing, data = df, clusters = mkid14) %>% tidy()
  
  bind_rows(mod1, mod2, mod3, mod4) %>%
    mutate(bw = bw,
           placebo = -10000) %>%
    filter(term == "opium_legal") %>%
    return()
}


bws <- c(2500, 5000, 7500, 10000)

placebo_positive <-
  lapply(bws, make_regs_bw_plus) %>%
  bind_rows()
  
placebo_negative <-
  lapply(bws, make_regs_bw_minus) %>%
  bind_rows()


placebo_plot <- 
  bind_rows(placebo_positive, placebo_negative) %>%
  mutate(bw_label = paste0((bw/1000), "km"),
         bw_label = factor(bw_label, levels = c("2.5km", "5km", "7.5km", "10km")),
         bw_label = fct_rev(bw_label)) %>%
  mutate(outcome_label = case_when(outcome == "outcome1" ~ "Place of worship",
                                   outcome == "outcome2" ~ "Live in village",
                                   outcome == "outcome3" ~ "Live in neighborhood",
                                   outcome == "outcome4" ~ "Rent room",
                                   TRUE ~ NA_character_)) %>%
  mutate(placebo_label = case_when(placebo == 10000 ~ "Placebo Boundary: +10km",
                                   placebo == -10000 ~ "Placebo Boundary: -10km")) %>%
  #mutate(sig_flag = (p.value < 0.1)*1) %>%
  ggplot(aes(x=estimate, y = outcome_label, group = bw_label, shape = bw_label, color = bw_label)) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbarh(aes(xmin=estimate-1.67*std.error, xmax = estimate+1.67*std.error), position = position_dodge(width = 0.4), height = 0) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  xlab("Effect of Opium Concession System") +
  #scale_color_manual(values = c("lightgrey", "black")) +
  facet_wrap(placebo_label~.) +
  guides(color = FALSE) +
  theme_bw() +
  scale_color_grey() +
  theme(axis.ticks.y.left = element_blank(),
        strip.background = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title.y = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.title = element_blank(),
        axis.line.y.left = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))

ggsave("./_4_outputs/figures/figure_6.tiff", plot = placebo_plot,  width = 5, height = 3, dpi = 300, compression = "lzw")
